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Abstract 

Hybrid shell elements have long been regarded with reserve by the commercial finite 
element developers despite the high degree of reliability and accuracy associated with 
such formulations. The fundamental reason is the inherent higher computational cost 
of the hybrid approach as compared to the displacement-based formulations. How- 
ever, a noteworthy factor in favor of hybrid elements is that numerical integration 
to generate element matrices can entirely be avoided by the use of symbolic inte- 
gration. In this paper, the use of the symbolic computational approach is presented 
for an assumed-stress hybrid shell element with drilling degrees of freedom and the 
significant time savings achieved is demonstrated through an example. 
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Introduction 


Evaluation of the finite element quantities (i.e., stiffness and mass matrices, load 
vector, etc.) is usually performed using numerical integration with an appropriate 
Gaussian quadrature rule (e.g., see [1]). However, two significant issues arise. First, 
if the integrand is a polynomial then an appropriate order Gaussian quadrature rule 
will result in exact integration. Second, as the complexity of element formulation 
increases so does the computational effort to evaluate element matrices and vectors. 
In the standard displacement formulation for parallelogram-shaped elements, the or- 
der of the Gaussian quadrature rule can be specified a priori. However for general 
quadrilateral-shaped (i.e., distorted) elements, the Jacobian is dependent on the el- 
ement natural coordinates, and hence, the integrand is a ratio of polynomials. As 
such, numerical integration with a given Gaussian quadrature rule results in an ap- 
proximation. 

In finite element formulations based on multi-field variational principles (mixed 
formulations), the element matrices could be generated explicitly such that element 
operations like numerical integration are entirely avoided. This is achieved by using 
symbolic computational techniques and symbolic manipulation software (e.g., see 
MACSYMA [2], MAPLE [3], MATHEMATICA [4]). 

The applicability of these computational techniques is by no means restricted to 
finite element analysis or computational mechanics alone (e.g., see Beltzer [5] and Tsai 
and Kikuchi [6]); however, as this paper deals with a successful implementation of a 
symbolic computational approach in a finite element code, in what follows we shall 
focus primarily on symbolic methods within the context of a finite element method. 
Symbolic and algebraic manipulation software systems for the evaluation of finite el- 
ement arrays has been used for approximately two decades beginning with the early 
work of Gunderson and Cetiner [7] and Luft et al. [8]. Andersen and Noor [9], Korn- 
coff [10], Korncoff and Fenves [11], Noor and Andersen [12, 13], and Tan et al. -[14] 
give an account of the development of symbolic manipulation in finite element code 
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generation. In each case the symbolic approach is noted for its robustness, accuracy, 
reliability, and for assisting in algebraic manipulations involved in research towards 
new finite element formulations. While most studies concentrated on displacement- 
based formulations, several researchers have examined mixed and hybrid formulations 
(e.g., Noor and Andersen [12], Saleeb and Chang [15], Chang et al. [16], and Tan 
et al. [14]). These studies focussed primarily on the symbolic computational aspects 
of the problem and present very limited data, if any, comparing the computational 
effort required for the symbolic and numerical approaches. Performance data com- 
paring other shell element implementations on a nearly level playing field in terms of 
computational effort to achieve a specific solution error have not been presented* 

The objective of this paper is to present comparison demonstrating the effective- 
ness of a symbolic computational approach over the numerical integration approach 
for a 4-node assumed- stress hybrid shell element with drilling degrees of freedom. 
In general, hybrid formulations involve more computations at the element level than 
displacement-based formulations and hence are often shunned by developers of com- 
mercial finite element software systems. However, hybrid formulations readily lend 
themselves to symbolic computational techniques even for distorted elements. The 
symbolic approach offers significant computational savings over the traditional nu- 
merical integration approach and thereby should make these elements attractive from 
both a computation and accuracy point of view. A brief description is first given for 
the hybrid element formulation and finite element approximations, followed by a de- 
scription of the use of the symbolic computational strategy. Finally, the effectiveness 
of the present strategy is demonstrated through a timing-study on a specific example. 

Hybrid Element Formulation 

The hybrid element formulation used in this paper is based on the Hellinger- 
Reissner multi-field variational principle and is given in detail in reference [17]. A 
brief outline of the formulation is given herein for completeness and to clarify the 
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notation used herein. The ensuing discussion pertains to a solid elastic continuum 
V, with boundary S. Let S„ be the section of the boundary where tractions are 
prescribed; let be the section of the boundary where displacements are prescribed. 
The Hellinger-Reissner functional can be written as 

n„* = -If {<t} t {d]{o} iv + f jv - / {u} T {(.} ds , (1) 

Z JV JV J Sff 

where [jD] is the compliance matrix, [£] is the linear differential operator on the 
displacements {u} to produce strains, and {<„} is the prescribed tractions on S„. 

The approximations for stresses and displacements can now be incorporated in 
the functional. The stress field is described in the interior of the element as 

M = [JW>, (2) 

and a compatible displacement field is described by 

M = [«]{,}, (3) 

where [P] and [ N ] are matrices of stress and displacement interpolation functions, 
respectively, and {/?} and {<?} are the unknown stress and nodal displacement param- 
eters, respectively. Substituting the stress and displacement approximations [equa- 
tions (2) and (3)] in the functional 11//;* [equation (1)] results in 

n„ B = w 

where 

m = J y [p) T mp]dv, 
m = J v [p] T [cm dv , 

{F} = / (5) 

Imposing stationary conditions on the functional with respect to the stress parameters 

{(}} gives 

{/?} = [H]-'[T}iq). ( 6 ) 


3 



( 7 ) 


Upon substituting equation (6) into equation (4), the functional reduces to 

n hr = — {9} T {7'’}> 

where [K] 9 the elemental stiffness matrix, is given by 

[I<] = [T] t [H)- 1 [T]. (8) 

Imposing the stationary conditions on the functional with respect to nodal dis- 
placement unknowns {q} results in a system of equations for the element of the form 

mm = ,< 9 ) 

Assembly of the element equations for a given finite element model and application 
of boundary conditions then defines the global system of linear algebraic equations 
to be solved. 

Displacement Field Approximations 

The generalized displacements and rotations represented in terms of interpolation 
functions and nodal values of the displacements and rotations are 

v°(t,ri) 

w°(Lv) 

Qx{Z,v) 
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and 


t + 1 ; i — 1,2, 3 


(13) 


r 

j = < 

1 ; i = 4 

In equations (10)— (13), t and tj denote the element natural coordinates, and ti and 
rji denote the values of t and rj at node i. Note that the inplane displacement field 
approximations include contributions from the drilling degrees of freedom in order to 
increase the level of approximation over the element without introducing additional 
nodes in the element. Further details on the displacement field and shape functions 
are given in references [17, 18). 

Stress Field Approximations 

The stress resultant field for the membrane part is assumed to be 

Ni = 0 1 + 0 4 t + PeV + PsV 2 , 

N v = 02 + Pit + 07V + 0st\ 

Nfr = 03 - 04V - 07t • O 4 ) 

The stress resultant field approximation for the bending part is assumed to be 

M(_ = 01 + 0Ai + 06V + 06V* i 

Mr, = 02 + 0st + 07V + Pot^ i 

M{„ = 03 + Piot + Pi iV + 2 ^ 12 ^ + 2 ^ 13 ^* 

The transverse shear stress resultants are obtained by relating the stress parame- 
ters of the transverse shear stress resultants to those of the bending stress resultants, 
as 

Q( = 04 + 011 + 013V> 

Qr, = 07 + 010 + Put- (1®) 

Note that while the membrane and bending stress resultant field approximations are 
uncoupled, the transverse shear stress resultant field is coupled to the bending stress 
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resultant field, and the equilibrium equations are satisfied a priori for regular (rect- 
angular) shapes. The stress field approximations are expressed in the element natural 
coordinate basis, and the stress interpolation function matrix [.P] is not uncoupled. 
This makes it virtually impossible to employ partitioning techniques in the evaluation 
of [i7] and [T] matrices a s was done in reference [14]. 

Symbolic Computation Approach 

The evaluation of elemental stiffness matrix [/f], equation (8), requires the eval- 
uation of [ H ] and [T] matrices [equation (5)]. The [H] matrix is rewritten in a 


convenient form as 


[H) = 

(17) 

where 

|tf«,->)] = det|J|. 

For this shell element, the matrix [/f] is a 22 x 22 matrix. 
Similarly the [T] matrix is expressed as 

(18) 

in = /'£[?({,>>)]#<(>) 

(19) 

where 

|f ({,»)) = [Pf |£][JV] det|J|. 

(20) 


For this shell element, the matrix [T] is a 22 x 24 matrix. 

In the symbolic evaluation of the [H] and [T] matrices, the integration is per- 
formed exactly and no approximation is introduced. Moreover the need to perform 
the integration for every element during an actual analysis is overcome - only the eval- 
uation of symbolically generated results is required. Note that in equation (20), the 
differential operator [£] acting on [N] implicitly introduces a det| J\ in the denomina- 
tor which cancels the det | J| appearing in the numerator. This is typical of the hybrid 
approach and makes it ideal for symbolic evaluation of elemental stiffness matrix. 
However in displacement-based elements, a det|J| will remain in the denominator 
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and hence will make direct symbolic integration very cumbersome or impractical, at 
least with the available software today, unless the element is restricted to a parallel- 
ogram shape with evenly spaced nodes and straight edges. One approach proposed 
by Andersen and Noor [9] to alleviate this inherent difficulty is to combine the use of 
symbolic and numerical integration methods. 

In the numerical integration version, a 3 x 3 Gaussian rule is used to integrate 
these matrices exactly. Symbolic evaluation of the matrix inverse of [f/] for this shell 
element leads to symbolic expressions that are too cumbersome for most FORTRAN 
compilers. This occurs even when the matrix inverse is represented by independent 
expressions for its determinant and cofactor matrix. Instead, the term [A] = [i/] -1 [T] 
is evaluated numerically by solving the following system of linear algebraic equations 
which has multiple right-hand-side vectors (i.e., columns of the matrix [T]), or 

[H)[A] = [T] (21) 

Hence the explicit matrix inversion is avoided and efficient, optimized equation solvers 
can be exploited. 

Results 

To demonstrate the effectiveness of the symbolic integration approach over the 
numerical integration approach, a timing study is presented. For convenience, the nu- 
merical integration version of the present element formulation is referred to as A4N1 
and the symbolic integration version as A4S1. These elements are implemented in 
the Computational MEchanics Testbed COMET using the generic element proces- 
sor [19]. As such, a common basis for the evaluation of different element formulations 
and implementations is readily and systematically accomplished. Common utilities 
for model generation, assembly and solution are available, thereby allowing element 
developers to focus only on element related issues. The other shell elements available 
in COMET include a 4-node incompatible assumed natural-coordinate strain element 
(4_ANS), a 4-node incompatible displacement-based element with drilling degrees of 
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freedom (4-STG), and a 4-node hybrid shell element without drilling degrees of free- 
dom (4-HYB). Note that the element routines for both 4.ANS and 4.STG have been 
highly optimized within COMET, whereas the element routines for both A4N1 and 
A4S1 are research-oriented implementations. With the exception of the element rou- 
tines, all other aspects of the analysis (e.g., data handling, element assembly, equation 
solving, etc.) are identical in that COMET provides a level playing field for assessing 
element performance. 

The emphasis here is to compare the use of symbolic and numerical integration 
procedures and to determine their effect on computation time for evaluating element 
matrices and vectors. As such, only a single example problem is considered. For this 
purpose, an isotropic square plate with all edges clamped and subjected to a central 
concentrated load is considered (see Figure 1 and reference [17]). Due to symmetry, 
finite element models of only one-quarter of the plate are used in the computations. 
The mesh is refined and the CPU time to form the stiffness matrix is calculated at 
each level of mesh refinement. 

The average computational effort to evaluate the element stiffness matrix as a 
function of the number of elements in the finite element model is shown in Figure 1. 
The CPU time per element converges to a constant value in all cases after approx- 
imately 30 elements. These results indicate that the 4.ANS and 4.STG elements 
require approximately the same computational effort and are the fastest of the el- 
ement implementations considered herein. Also note that the symbolic integration 
approach (A4S1) requires only half the computational effort of the numerical integra- 
tion approach (A4N1). In effect, the symbolic integration approach offers a gain of 
almost 50% computational time over the numerical integration approach. It is antici- 
pated that another 10-20% reduction in CPU time per element is achievable once the 
element routines are optimized. Currently A4S1 requires just under twice the CPU 
time per element as either 4_ANS or 4-STG to evaluate the linear stiffness matrix. 
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A performance comparison should also take into account the accuracy of the 
predicted solution by different elements. It is possible to define a weighted average 
performance measure that takes into account both the CPU time and the error in 
the solution. The real difficulty lies in assigning appropriate weights for CPU time 
and solution error, since these depend on the user requirements. Hence in the present 
study, a more general representation is adopted by comparing the CPU time with 
the solution error. Solution error is defined for this problem as the absolute value of 
•the ratio of the difference between the value of the center transverse deflection from 
classical plate theory and the value from the finite element analysis to the classical 
plate theory value, in percent. 

The variation of the solution error as number of elements increases is shown in 
Figure 2. It is evident from Figure 2 that the present element formulation imple- 
mented as either A4N1 or A4S1 offers the best solution in that it converges to the 
exact value faster than the other elements considered. Using the results given in 
Figures 1 and 2, the CPU time to form all the element stiffness matrices of the finite 
element mesh may be estimated. Table I gives an indication of the computational 
effort to achieve a specified solution error. The total CPU time to evaluate all the el- 
ement stiffness matrices (t K ) and the overall solution time (<°) for each element type 
considered are normalized with respect to the corresponding times obtained using the 
A4S1 element and are presented in the last two columns of Table I. For a specified 
solution error percentage, 4_ANS takes the least time. Among the other elements, 
A4S1 offers the best time. In fact, 4_STG and 4_HYB do not predict a solution error 
of less than 0.5 % even with 144 elements. Though 4_ANS uses a larger number of 
elements to obtain a specified solution error than A4S1, it offers a better CPU time 
since it requires only half the CPU time per element compared with A4S1. However, 
both 4.ANS and 4_STG have been highly optimized for performance in terms of CPU 
time in COMET, while the element routines for A4S1 axe still in a development form. 
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In general, as the robustness of an element formulation increases, the computa- 
tional effort required to evaluate element properties also usually increases. For these 
hybrid elements, symbolic computations significantly reduce the CPU time required 
to evaluate the element stiffness matrix (compare results for A4S1 and A4N1). Ad- 
ditional improvements in the computing times are possible by reviewing the symbol- 
ically generated FORTRAN code and making any necessary enhancements to reduce 
the CPU time per element. Furthermore, an optimization study for A4S1 similar to 
4_ANS and 4_STG is expected to yield additional computational time savings. With 
these enhancements, it is anticipated that the A4S1 element will exhibit total CPU 
times comparable to the 4_ANS element and perhaps better. 

Furthermore, it is worth noting that both the 4_ANS and 4_STG elements are 
more sensitive to mesh distortion than the present A4S1 hybrid element as shown in 
reference [17]. While the performance data presented herein is for a uniform mesh 
of square elements, the computational effort to achieve a specific solution error with 
a mesh of distorted elements will increase significantly for the 4-ANS and 4-STG 
element cases. In these cases, the number of these elements required to achieve a given 
accuracy will be much larger than that required for the A4S1 element. Therefore, it 
is anticipated that for distorted meshes, the A4S1 hybrid element will offer improved 
performance over the 4_ANS element. 

Conclusions 

The use of symbolic computational approach in a 4-node hybrid shell element 
is presented and is shown to be almost twice as fast as the numerical version of the 
same element. This gain of computational time is significant as it makes the hybrid 
element, hitherto known to be computationally costly despite being more accurate, 
compete with efficient 4-node displacement-based elements. 

On the symbolic approach front, the MAPLE output expressions are rather un- 
wieldy; however, the expressions are automatically converted to FORTRAN state- 
ments and hence could be directly incorporated in to an existing finite element code. 
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Table I. Computational effort required to achieve a specified solution error. 


Solution 

Element 

Number of 



Error 

Name 

elements 

i K li K A 451 

<5 % 

A4S1 

9 

1.000 

1.000 


A4N1 

9 

1.450 

1.044 


4_ANS 

16 

0.763 

0.835 


4.STG 

16 

0.919 

1.098 


4.HYB 

16 

1.850 

1.146 

<2 % 

A4S1 

16 

1.000 

1.000 


A4N1 

16 

1.650 

1.114 


4_ANS 

25 

0.793 

0.883 


4_STG 

36 

1.147 

1.375 


4.HYB 

49 

3.622 

1.738 

<1% 

A4S1 

36 

1.000 

1.000 


A4N1 

36 

1.783 

1.198 


4_ANS 

49 

0.758 

0.891 


4.STG 

100 

1.440 

1.978 


4.HYB 

121 

4.486 

2.499 

< 0.5 % 

A4S1 

64 

1.000 

1.000 


A4N1 

64 

1.809 

1.269 


4.ANS 

81 

0.706 

0.884 


4.STG 

- 

- 

- 


4.HYB 

- 

- 

— 


t K : Total CPU time to evaluate all the element stiffness matrices. 
t 1 ® : Overall solution time (total time from start to finish). 

t^i : Total CPU time to evaluate all the element stiffness matrices using A4S1. 
t° 4S1 : Overall solution time using A4S1. 




